cL2p_ls_tsvm <- function(X, y, C1 = 1, C2 = C1, 
                         kernel = c("linear", "rbf", "poly"),
                         gamma = 1 / ncol(X), degree = 3, coef0 = 0,
                         p = 1, epsilon1 = 1, epsilon2 = epsilon1,
                         eps.irls = 1e-2, irls.steps = 10,
                         fit_intercept = TRUE, randx = 1) {
  cL2p_ls_tsvm_dual_solver <- function(KernelX, idx, C1, C2, 
                                       p, epsilon1, epsilon2,
                                       eps.irls, irls.steps) {
    update_weight <- function(X, w, p, epsilon) {
      f <- X %*% w
      f2 <- f^2
      fL2p <- f2^(p/2)
      f2[(f2 < 1e-7)] <- 1e-7
      weight_elem <- (p/2)*f2^((p - 2)/2)
      weight_elem[fL2p >= epsilon] <- 0
      return(as.vector(weight_elem))
    }
    
    X1 <- KernelX[idx, ]
    X2 <- KernelX[-idx, ]
    
    xn <- nrow(KernelX)
    xp <- ncol(KernelX)
    X1n <- nrow(X1)
    X2n <- xn - X1n
    
    Fvec <- rep(1, X1n) 
    Dvec <- rep(1, X2n) 
  
    epsI <- diag(1e-7, xp)
    
    e1 <- matrix(1, X1n)
    e2 <- matrix(1, X2n)
    
    u01 <- matrix(0, xp)
    u02 <- matrix(0, xp)
    
    GramX2 <- t(X2) %*% X2
    GramX1 <- t(X1) %*% X1
    
    X1T_e1 <- t(X1) %*% e1
    X2T_e2 <- t(X2) %*% e2

    for (i in 1:irls.steps) {
      u11 <- -cholsolve(t(X1*(Fvec/C1)) %*% X1 + GramX2 + epsI, X2T_e2)
      Fvec <- update_weight(X1, u11, p, epsilon1)
      if (norm(u11 - u01, type = "2") < eps.irls) {
        break
      } else {
        u01 <- u11
      }
    }
    
    for (i in 1:irls.steps) {
      u12 <- cholsolve(t(X2*(Dvec/C2)) %*% X2 + GramX1 + epsI, X1T_e1)
      Dvec <- update_weight(X2, u12, p, epsilon2)
      if (norm(u12 - u02, type = "2") < eps.irls) {
        break
      } else {
        u02 <- u12
      }
    }
    
    BaseDualCL2PLSTSVMClassifier <- list("coef1" = as.matrix(u11),
                                         "coef2" = as.matrix(u12))
    return(BaseDualCL2PLSTSVMClassifier)
  }
  X <- as.matrix(X)
  y <- as.matrix(y)
  class_set <- sort(unique(y))
  idx <- which(y == class_set[1])
  y[idx] <- -1
  y[-idx] <- 1
  y <- as.matrix(as.numeric(y))
  if (length(class_set) > 2) {
    stop("The number of class should less 2!")
  }
  kernel <- match.arg(kernel)
  solver <- "dual"
  kso <- SupportVectorLab:::kernel_select_option(X, kernel, "primal", randx,
                                                 gamma, degree, coef0)
  KernelX <- kso$KernelX
  Kw <- kso$KernelX[kso$sample_idx, ]
  X <- kso$X
  if (fit_intercept == TRUE) {
    KernelX <- cbind(KernelX, 1)
  }
  if (solver == "dual") {
    solver.res <- cL2p_ls_tsvm_dual_solver(KernelX, idx, C1, C2,
                                           p, epsilon1, epsilon2,
                                           eps.irls, irls.steps)
  }
  TSVMClassifier <- list("X" = X, "y" = y, "class_set" = class_set,
                         "C1" = C1, "C2" = C2, "kernel" = kernel,
                         "gamma" = gamma, "degree" = degree, "coef0" = coef0,
                         "solver" = solver, "coef1" = solver.res$coef1,
                         "coef2" = solver.res$coef2,
                         "fit_intercept" = fit_intercept,
                         "Kw" = Kw)
  class(TSVMClassifier) <- "TSVMClassifier"
  return(TSVMClassifier)
}